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Unified Framework and Algorithm for Quantized Compressed Sensing 

Zai Yang*, Lihua Xie*, and Cishen Zhang^^ 



Abstract 



Compressed sensing (CS) studies the recovery of high dimensional signals from their low dimensional linear 
^D measurements under a sparsity prior. This paper is focused on the CS problem with quantized measurements. 
There have been research results dealing with different scenarios including a single/multiple bits per measurement, 
noiseless/noisy environment, and an unsaturated/saturated quantizer. While the existing methods are only for one 
or more specific cases, this paper presents a framework to unify all the above mentioned scenarios of the quantized 
<^ CS problem. Under the unified framework, a variational Bayesian inference based algorithm is proposed which is 
pg applicable to the signal recovery of different application cases. Numerical simulations are carried out to illustrate 
PsJ the improved signal recovery accuracy of the unified algorithm in comparison with state-of-the-art methods for 
both multiple and single bit CS problems. 



1 Introduction 



The recently developed compressed sensing (CS) theory and methods |[T]|2| can achieve acquisition of information 
contained within a huge volume of data using only a small amount of measurement samples. Different from the 
classical Shannon-Nyquist sampling theorem which requires that the sampling frequency be twice as many as the 
f^ bandwidth of a signal in order to reconstruct its complete information, the CS theory assesses the success of signal 

OO recovery with the sparsity. A signal x G M^ of length N is called K-sparse in a basis ^ G M^^^ if all but at 

^ most a number of K <^ N entries of its coefficient 6 € M^ are zeros with x = ^6. Without loss of generality 

^ we assume that * is an identity matrix, i.e., x is sparse in the canonical basis, since for a general basis * it can be 

^ absorbed into the following introduced sensing matrix A. Rather than observing directly the original sparse signal 

^—1 X, a number of M, K < M <^ N, linear measurements are acquired in CS as 

r^ y = Ax + n (1) 

X 

5^ where y G M^^ is the measurement vector, A G M^^^^ denotes the sensing/measurement matrix and n G M^^ is 

the measurement noise vector. Though the recovery of x from y is generally ill-posed (less linear equations than the 
unknown variables), it is shown in [3] that a sparse signal x can be stably recovered under mild conditions on A in 
the sense that the recovery error grows linearly with the noise level. To do this, the basis pursuit denoising (BPDN) 
problem 

min||i||j^, subject to ||y — A;r||2 < e (2) 

is solved where e > ||n||2 denotes the noise level. The recovery is exact in the noise free case. A similar result holds 
for compressible signals that are not exactly sparse. 

Sparse Bayesian learning (SBL) ||4|[5| was derived from the research area of machine learning and has become 
a popular method for sparse signal recovery in CS. In SBL, the sparse signal recovery problem is formulated from a 
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Bayesian perspective while the sparsity information is exploited by assuming a sparse prior for the sparse signal of 
interest. As an example, a Laplace prior corresponds to the ii norm widely used in optimization approaches ||5j[6|. 
Since the exact Bayesian inference is typically intractable, approximation approaches to Bayesian inference have 
been adopted including evidence procedure [7|, e.g., in |j6|, and variance message passing (VMP) |81, e.g., in Q. 
One merit of Bayesian CS is the flexibility of modeling sparse signals that can not only promote the sparsity of 
its solution, e.g., in \^, but also exploit additionally known structure of the sparse signal, e.g., in [10]. Since the 
Bayesian inference is a probabilistic method and based on heuristics to some extent, one shortcoming of Bayesian 
approaches is that there have been fewer results on their signal recovery accuracy in comparison with deterministic 
approaches, e.g., BPDN. 

The conventional CS framework is mainly focused on the sparse signal recovery from the real-valued measure- 
ment y that has infinite bit precision. The required number of measurements M is mainly studied for guaranteed 



signal recovery accuracy pT| - [T4| . Since quantization is necessary for practical considerations, e.g., data storage and 
transmission, we study the sparse signal recovery from quantized measurements in this paper. During the quanti- 
zation process, each continuous-valued measurement is quantized into some value in a finite set. A new challenge 
is thus the existence of quantization errors. The noise free case with a uniform unsaturated quantizer is studied 
in |[3 15 16 1. A solver with quantization consistency is recommended in 1 15 1 that corresponds to replacing the I2 



norm in BPDN by the l^o norm. A BPDN solver is used in [3] to handle the quantization errors as additive white 
Gaussian noises because of the stable signal recovery performance of BPDN. A family of solvers, named as basis 



pursuit dequantizer of moment p (BPDQp), that includes BPDN and that in 1 15 1 as special cases is studied in 1 16 1 
where the (.2 norm in BPDN is replaced by an ip norm with 2 < p < +00. By characterizing the quantization errors 
as independent random variables uniformly distributed in a common interval, it is shown in |16| that the optimal 
signal recovery accuracy is obtained at some finite p > 2. But unfortunately, the optimal p cannot be explicitly given 
in practice. Note that the common uniform distribution assumption is crucial to obtain the results in p6| . As a result, 
it is unclear whether the results in [ 16] can be extended to a general quantizer case where such an assumption fails. 
It is obvious that both BPDN and BPDQp are inappropriate in the case of a saturated quantizer since data saturation 
may lead to large or even unbounded quantization errors that deteriorate their performance. To deal with the data 
saturation, Laska et al. ||T7| propose two modified versions of BPDN, including saturation rejection and consistency 
approaches, denoted by BPDN-SR and BPDN-SC respectively. The saturated measurements are simply rejected in 
BPDN-SR while they are incorporated in the signal recovery process in BPDN-SC. While quantization errors and 
measurement noises are coupled in most existing methods (some methods, e.g., BPDQp, consider only the noise free 



case to avoid such a problem), e.g., in 1 17 1, they are separately studied by Zymnis et al. [18] where the authors seek 
to find a signal estimate that maximizes the likelihood of the quantized measurements while ii norm is exploited 
to promote the signal sparsity. The resulting algorithm is quoted as £1 -regularized maximum likelihood (LIRML) 
algorithm. 

An extreme case of quantized CS is so-called 1-bit CS where each quantized measurement keeps only the sign 



information of the real- valued measurement and thus uses just one bit. The 1-bit CS framework is proposed in j 19 1 
and has attracted many research interests because it possesses many merits. For example, a 1-bit quantizer is a simple 
comparator that tests whether the measurement is above or below zero, leading to an easy implementation and a fast 
quantization process. A measurement noise can be neglected in 1-bit CS as long as it does not change the sign of 



the measurement. It is shown in [20 1 that to acquire just one bit for each measurement is optimal in the presence of 
heavy noises. The 1-bit case is quite different from the multi-bit case since all measurements are saturated in 1-bit 
CS and the signal scaling information is lost. A common approach to the signal scaling problem is to impose that 
the signal to be recovered has a fixed unit norm and then search for the signal on the unit hyper-sphere rather than in 
the whole space. Such a constraint is nonconvex and brings new challenges to algorithm design. Existing algorithms 
based on this constraint include renormalized fixed point iteration (RFPI) [19], matching sign pursuit (MSP) | |2T) , 
restricted-step shrinkage (RSS) [22] and binary iterative hard thresholding (BIHT) [|23j|. Convex formulations of 



the 1-bit CS problem have been recently proposed by Plan and Vershynin ||24 25 1. They show in p4J that a linear 



program can decode the noiseless case with guaranteed signal recovery accuracy under similar mild conditions as in 



conventional CS. In flSl they introduce a seemingly unrelated convex program for the noisy case and show similar 
results. It is noted that both the BIHT and the convex problem in | ,25J that deal with the noisy 1-bit CS problem 
require the signal sparsity information (BIHT needs to known ||a;||Q and CVXP requires a proper upper bound for 
IItII t 

In this paper, we introduce a unified framework for quantized CS that works for both multi- and 1-bit CS cases. 
The new framework deals with quantization errors and measurement noises separately, allows data saturation in 
multi-bit CS, and does not need the signal sparsity information. Based on Bayesian inference, we show that many 
existing formulations for quantized CS are special cases of the new framework or related to it after minor mod- 
ifications. Based on the new problem formulation, we propose an algorithm within the Bayesian CS framework 
where a three-layer hierarchical prior introduced in [9] is adopted as the sparse signal prior and Bayesian inference 
is carried out using VMP. The performance of the proposed algorithm is studied by extensive numerical simulations 
in various scenarios. It is shown that the new algorithm improves the signal recovery accuracy in comparison with 
state-of-the-art methods in both multi- and 1-bit CS. 

Notations used in this paper are as follows. Bold-face letters are reserved for vectors and matrices. For ease of 
exposition, we do not distinguish a random variable from its numerical value. Xi is the ith entry of a vector x. xx 
denotes a truncated vector of x with entry indices in a set I. Ai is the ith column of a matrix A. \\x\\q counts the 
number of nonzero entries of a vector x. \\x\\ = (^^ \xif) denotes the Ip norm of a vector x with 1 <p < +00. 
{g {x)) /^\ denotes the expectation of a function g {x) with respect to a random variable x whose probability density 
function is p(x). denotes Hardamard (elementwise) product. )^ and ^ denote > and < respectively with an 
elementwise operation. 

The rest of the paper is organized as follows. Section [2] introduces the unified framework for quantized CS 
and studies its relations with existing formulations. Section[3]introduces the proposed Q-VMP algorithm. Section|4] 
presents numerical simulations to illustrate the improved signal recovery accuracy of the proposed Q-VMP algorithm 
in comparison with existing ones. Section [5] concludes the paper and discusses some future works. 

2 A Unified Framework for Quantized CS 

2.1 A Unified Observation Model 

In quantized CS, the observed samples are noisy linear measurements of the original signal after quantization: 

z = Q{y), y = Ax + n (3) 

where x is the signal of interest, A is the sensing matrix, n is the measurement noise vector, y is the pre-quantized 

noisy measurement vector, Q denotes a quantizer and z is the observation. A (regular) quantizer Q {v) for a scalar- 

u G M is defined as 

vq, if w G (no,ni), 

Q{v)=\ ""][ if;^eK,n2), ^4^ 

' ' ' 1 * * * ; 

^ VL-l, if u G [ul-i,ul), 

where L denotes the quantization level and typically satisfies L = 2^ with B denoting the bit depth (bits per 
quantized measurement), uq < ni < • • • < ul, and Vi G [uj, Uj+i) for i = 0, • • • , L — 1. The quantizer Q (v) is 
called unsaturated if {uq, ul) is a finite interval, or saturated otherwise. For a vector v, Q (v) operates elementwise. 
Multi -bit CS refers to the case B > 2 while 1-bit CS corresponds to B = 1. 

2.1.1 Multi-bit CS 

We consider first a multi-bit quantizer where B > 2. Denote Vy the domain of y. Then we have 

Vy = Q-'{z):={yGR^'\Qiy) = z}. (5) 



We introduce an auxiliary variable e denoting the quantization error: 

e = z-y 

with its domain 

Ve = z-Vy:={z- y\y G Vy] . (6) 

Note that Pg is unbounded as data saturation occurs. 

2.1.2 1-bit CS 

In the case of 1-bit quantizer we set uq = —oo,ui = and U2 = +oo. The sign information of y is preserved in 
the quantized measurement z. But the scaling information of y and that of x is lost. Without loss of generality, we 
let the 1-bit quantizer 

Q (v) = ?sgn (v) 

for a scalar u G M with ij — )• 0+ (<j is an arbitrarily small positive number) and sgn (•) being the sign function. For 
convenience, we set sgn (0) = 1 (the choice is arbitrary and can be replaced by sgn (0) = —1). Then we have 
2; — ;■ 0. To solve the signal scaling problem we impose a constraint that y has fixed unit norm, i.e., 

Ili/IL = 1 (7) 

with s > 1. Different from the multi-bit quantizer case we have in such a case that 

Vy = {yG M^lsgn (y) = sgn (z) , \\y\l = l} , (8) 

Ve = {ee M^-^lsgn (e) = -sgn (z) , ||e||, = l} . (9) 

As a result, a unified model for both multi- and 1-bit CS can be written into 

z = Ax + e + n, e G Pg, (10) 

which is the observation model to be used in this paper to recover x. 



Remark 1 In the observation model (10), y is written into two variables z and e that are dependent. As an 
example, the statistical information of e can be obtained given that of x and n for a fixed quantizer Q. For the 
sake of computational ease, such dependence will not be exploited in the following signal estimation process, i.e., 
we assume that there is no a priori information ofe but its domain T)^,. 

2.2 Bayesian Formulation of the Quantized CS Problem 

In this subsection we formulate the quantized CS problem from a Bayesian perspective based on the observation 
model in ( 10 1. According to Remark[l]we treat e as a variable independent of x and n. The joint probability density 



function (PDF) p {z, x, e) is decomposed as 

p{z,x,e) =p{z\x,e)p{x)p{e). 
We define the three distributions on the right hand side as follows. 

2.2.1 Noise model 

Under an assumption of white Gaussian measurement noise, i.e., n ~ A^ (O, o-"^!) where a"^ is the noise variance 
and / denotes an identity matrix of proper dimension, we have 

p(z\x,e;a'^) = M (z\Ax + e,a'^l) . (11) 



2.2.2 Sparse signal model 

A sparse prior is needed for the sparse signal x of interest. Here we do not give an explicit distribution to the sparse 
signal X but denote p [x) its PDF. Then we let / [x) = —Ci \ogp{x) + C2 where Ci and C2 are proper constants. 
The only thing that we assume for p {x) is that it favors entries of x being zeros. As an example, a commonly used 
sparse prior for ic is a Laplace prior: p {x) = \^ exp {— A ||a3||x} with A being a positive constant. In such a case, 
we have / (a;) = \\x\\-y that has been extensively studied in deterministic optimization methods. 

2.2.3 Quantization error model 

We assume a uniform, noninformative prior for e: 

e^U {Ve) (12) 

since the only information of e that we use is e G "Dg. 

To obtain a maximum a posteriori (MAP) estimator of x requires to integrate out e from p {z, x, e) that is 
computational intractable. We propose to estimate x and e simultaneously using their joint MAP estimator: 

{x,e.} = arginaxlogp (a;, e\z) 

= argmaxlogp(2;,a;, e) 

x,e 

= argniaxlog{p (z|a;, e)p (a;)p(e)} ^ ' 



x,e£Ve 1 2a'^ 



arg min -^ / [x) + -^^-^ \\z — e — Ax\\2 



An equivalent form of the problem in ( 13 1 is 



niin/(x), subject to <j !! II2 ' ^^^-^ 



where e is a proper scalar that controls the noise energy. The first constraint in ( 14 1 is to ensure the data consistency 
against the measurement noise. In multi-bit CS, the second one concerns data consistency due to quantization. In 
1-bit CS, an additional signal scaling constraint is included in the second constraint that prevents an optimal solution 
for X from 0. Before proceeding to our algorithm within the framework of Bayesian CS, we study in the next 
subsection relations of the proposed Bayesian framework with existing methods. 

2.3 Relations with Existing Methods in Quantized CS 

We first note that problem ( [14] ) is equivalent to the problem 

min/(i), subject to \ 1^ ~ ^Il2 - ^' ^i^^ 

In the following we show that many existing problem formulations of quantized CS are special cases of or related to 



( 15 1, with the focus on the 1-bit case. 



2.3.1 Multi-bit CS 

We consider the case of £1 optimization where / {x) = \\x\\-^. In the noise free case where e = 0, the problem in 



( [15] ) can be written into 

min ll^ll^ , subject to Ax G Vy, (16) 



which has been studied in 1 26 1 . Further, by assuming that Q is a uniform unsaturated quantizer the above problem 

can be written into 

r 
min ^1 , subject to \\z — Ax\\^ < - (17) 

~ \\ \\ i_ J J II II CXJ — o 



that is studied in |T5}T6| with r denoting the quantization bin width. While existing methods that take into account 



measurement noise typically mix up the quantization error and the noise, e.g., in p7| , problem ( (T5| ) extends existing 
noise free formulations to the noisy case by deaUng with the two uncertainties separately. 



Remark 2 Under the assumption that all quantization errors are independent and uniformly distributed in a com- 
mon interval [—5,5], it is shown in /|i6|/ that the i-00 norm in problem {17) is not the best choice for the signal 



recovery. But it is unclear whether the result in [16'] can be extended to the case of a general quantizer where the 
above assumption fails. It is noted that our problem formulation does not require this assumption and applies to an 
arbitrary quantizer That is, by losing some optimality, we have obtained the universality. 



2.3.2 1-bit CS 



In 1-bit CS problem ( 15 1 becomes 



\y-Ax\\^ < e, 
min / {x) , subject to ^ sgn {y) = sgn (z) 

|j/L = i- 



(18) 



x,y 



In the noise free case, it can be written into 



min / (x) , subject to 



sgn (Ax) = sgn (z) , 



|A;r||, = l. 



(19) 



This problem with the settings f (x) = \\x\\-^ and s = 1 can be shown to be convex and has been studied in |24 1. So 



by problem ( (T8| ) we extend ( (T9| ) to the noisy case while the authors of p4| state in p5| that "it was unclear how to 
modify the above convex program to account for possible noise." 



Since the third constraint in ( 18 1 serves to prevent the optimal solution for x from 0, it is natural to replace it by 



I a; II 2 = 1, leading to the problem 



\y-Ax\\2 < e, 
mill / (x) , subject to ^ sgn (y) = sgn (z) 

\x\\n = 1. 



(20) 



x,y 



In the noise free case it becomes 



min / (x) , subject to 



sgn (Ax) = sgn (z) , 



X 



1, 



(21) 



which with / (x) = \\x\\-^ is the earliest formulation of the 1-bit CS problem introduced in p9J and solved re- 



spectively using RFPI in 1 19] and RSS in [22]. Note that for any optimal solution {x*,y*) to (20l, {x*,y') with 
y' = sgn (z) (sgn (z) Ax*)^ is also its optimal solution by 

\\y' - Ax*\\^ < \\y* - Aa;*|L < e, 

where {v), = max{w,0} for a scalar v and operates elementwise for a vector. So we can set y = sgn(z) 



(sgn (z) Ax)_^_ in ( 20 1, leading to that problem ( 20 1 is equivalent to the problem 

"(sgn(z)0 Ai)_|| <e, 



min / {x) , subject to 



I*ll2 = 1 



(22) 



where {v)_ = min {v, 0} for a scalar v and operates elementwise for a vector. 

The noise information is used in above formulations that account for the measurement noise. If it is unknown 



but the signal sparsity information is known, ( 22 1 can be written into 



min ||(sgn (z) Ai)_|| , subject to 



/ (x) < S, 

11-7*11 — 1 



(23) 



X 



where the constant Prefers to the sparsity information. BIHT-^2 is proposed in |i23|l to solve (23 1 with / (x) — \\^\\q. 
Moreover, the objective function in ( 23 1 has been also studied in | [T9|21|27 1 . By replacing the £2 norm in the objective 
function in (23 1 by £1 norm (that corresponds to assuming a Laplace distributed noise rather than Gaussian) it gives: 



mm 

X 



|(sgn (z) Aa;)_|| , subject to 



f{x)<S, 

Wnr-W — 1 



(24) 



which with / (x) = \\x\\q is solved using BIHT in |23|. Problem (24l searches for a solution in a sparse region 



(defined by the feasible domain) that maximizes the data consistency by minimizing the sum of the negative part of 
sgn (z) Ax. A related (not equivalent) problem is thus to maximize the sum of the positive part of sgn (z) Ax, 
i.e., to maximize ||(sgn (z) Aa;)^|| subject to the same constraints. By combining the two objective functions, 
i.e., to maximize || (sgn (z) Ax] 
IItII '^ 



+ 1I1 



(sgn (z) Ai) I = sgn (z) Aa;, it gives the problem (with / (a;) 



whose convex relaxation given by 



max sgn"^ (z) Ax, subject to 



max sgn {z)Ax, subject to 



\x\ 
\x\ 



\x\ 
\x\ 



<S, 



<S, 
< 1 



(25) 



(26) 



has been studied in |25 1 and shown to give guaranteed signal recovery accuracy under mild conditions. We note that 



it is easy to show that problem (25 I has the same guarantee as (26 1 following the analysis in |25 



Remark 3 



(1) Problem {18) with f (x) = \\x\\^ and s = 1, and problem (26) are two parallel convex formulations of noisy 



1-bit CS. Problem ( |i8| ) requires the noise information but not the signal sparsity level while it is an opposite 
situation for problem ( 26 ). Since guaranteed signal recovery performance has been proven for ( 26 ), we believe 



it is true that a similar result holds for (18), which we pose as an open problem and is beyond the scope of this 
paper. 



(2) By (22 ) we see that the effective noise is {sgn (z) Ax)_ in 1-bit CS where by "effective noise" we refer to 



a noise that has the minimum energy and leads to the same measurement. So, the effective noise is generally 
a sparse noise whose sparsity level depends on the sign flip rate of the measurement. It explains why the noise 



can be modeled by a Laplace distribution as in ( 24 ). Moreover, it explains why BIHT outperforms BIHT-£2 in 



the high SNR region, as reported in l\23'j , where the effective noise is very sparse. 

(3) It can be shown that the energy of the effective noise {sgn (z) Ax)_ is much smaller than that of the true 
noise y — Ax. From this point of view, we may say that 1-bit CS is robust to the measurement noise. 



Figure 1: Directed graphical model that encodes the joint PDF in (30 1 of the Bayesian model. Nodes denoted with 
circles correspond to random variables, while nodes denoted with squares correspond to parameters of the model. 
Doubly circled z is the observation while single circled nodes represent hidden variables. 



3 Q-VMP: Variational Message Passing for Quantized CS 

3.1 Model Selection 

We assume that the noise variance a^ is known. Though it can be estimated by assuming an inverse Gamma prior for 
it in the case where it is unknown as in |'6','28|, its estimate is inaccurate due to an "identifiability issue" as addressed 
in 129). For the sparse signal x, we adopt a three-layer, Gaussian-Gamma-Gamma, hierarchical prior introduced 



m 



0: 



p(.;e,c,<i) = //M-HpH.;e)M.;c,^)<i«d, 



where 



'p{x\cx) = A/'(a;|0, A) , 



TV 



p{oL\ri\e) = J|r(ai|e,??) 



4 = 1 



p{iT,c,d) = T{r]\c,d) 



(27) 
(28) 
(29) 



with A = diag (a) and constants e, c, d. For a Gamma distributed variable n ~ F (c, d), its PDF is F {u\c, d) = 
"^ exp (—du) with F (c) being the Gamma function. By |9] the constants e, c, d satisfy that < e < 1, 



r(c) 

c, a > 0. In this paper, we adopt c = 1, d = to make the prior for r/ in (29) noninformative (flat on M_|_). Further, 

we choose e = since a smaller e leads to a sparser prior and an estimator that approximates a hard-thresholding 

rule according to |[9| . Readers are referred to Q for more properties of the Gaussian-Gamma-Gamma prior and its 

relations with other sparse estimation techniques. 

In 1-bit CS, we let y have unit £2 norm in (jvl), leading to that ||e||2 = 1 in (J9]l. As aresult, we have the joint PDF 

of the observation model (fTOl): 



p{z,x,e,a,r]) = p{z\x,e)p{x\a)p{a\rj)p{r])p{e) 



(30) 



with the distributions on the right hand side as defined respectively by ([TTjl, ( [27] ), ( [28| ), ( [29] ) and ( [12] ). A directed 
graphical model that encodes the factorization of the joint PDF in ( 30 ) is shown in Fig. [T] 



3.2 Q-VMP Algorithm 

It is known that Bayesian inference is based on the posterior distribution p {x, e, a, r]\z) = p (z, x, e, a, t]) jp (2;). 
However, such an exact posterior distribution is intractable since p{z) = f ■■■ f p{z,x,e,a,r]) dxde da drj can- 
not be expressed explicitly. 

A variational inference approach |30 31) is adopted in this paper. Denote V = {ic,e,Q,T?} the set of all 
unknown variables to be estimated. The goal in variational inference is to find a tractable distribution q {V) that 
closely approximates the true posterior distribution^ (Viz). To do this, some family of distributions that has enough 
flexibility is firstly chosen to represent q{V). Then the task is to find a member of the family that minimizes the 



Kullback-Leibler (KL) divergence between the true posterior p{V\z) and the variational approximation q{V). A 
commonly used variational distribution q (V) is such that disjoint groups of variables are independent, i.e., q (V) 
has a factorized form q {V) = q (x) q (e) q (a:) q {rj). Variational message passing (VMP) is proposed in |[8| for the 
variational inference using a message passing procedure on a graphical model. In VMP, the variational distributions 
q {x), q{e),q {ex), q {rj) are iteratively updated to monotonically decrease the KL divergence and thus has guaranteed 
convergence. Readers are referred to L8J for more details of VMP. The updates of q {x), q (a), q {if) are similar to 
those in |[9l because of the similarity between quantized and conventional CS. q (e) is given complete flexibility in 
multi-bit CS as q {x), q {ex) and q {rj). We constrain q (e) in 1-bit CS such that 



'{e) = 6{e-e^ 



(31) 



pM 



due to a computational issue as to be discussed in Remark Wl where 6 {■) is a delta function and e^ € M^" is to 
be estimated. Note that ( [3T] ) is equivalent to the complete flexibihty in the noise free case as to be illustrated in 
Subsection [331 



3.2.1 Updates of q {x), q {ex) and q {rj) 
According to |[8J we have that 

q {x) oc exp | (Inp {z\x, e))^(g) (Inp {x\ex)) ^^^^^j 
fx exp I -- (a; - /i)^ S"^ {x- fj,)\ , 
and thus q{x) is a Gaussian distribution AA (a;|/x, S) with the mean /x and co variance Xl: 



/g(e) 



'<?{«) 



(32) 
(33) 



For ex we have 



q {ex) oc exp | {Inp {x\ex))^^^^ {Inp (« |^))g(^) } 



N 



OC JJ On ^ exp 

n=l 



'2«n^(4>g(^)-«n(??)g(^) 



where (x^) , ^ = /x,^ + S„„. The expression on the right hand side is the product of generalized inverse Gaussian 



(GIG) PDFs and thus we have for any z £ R plj: 



{xl),,^-,V^.+i^ ' 



\. *^fr 



2 {'n)q{r]) {^n)q(x) 



g{a) 



2(7?) 



q(v) 



^e-l 



Y 2 (??)q(.^) {Xn)q{x) 



(34) 



where K.,^ {■) is the modified Bessel function of the second kind and order i/ G M. The case of i = — 1 in (34i gives 
the evaluation of (A^^) , s used in (33 1, and the case of i = 1 gives the calculation of {an)q(ct) ^^^^ ^^ ^ l^^er 

expression in (35 1. The update of q (r/) can be shown to be q{r]) = T (ri\N€ + c, X]„=i {'^n)q(a) + d). The first 



moment of r] used in (p4b is given as 



iv) 



Ne + c 



<i{v) 



^Af 



En=l («n)<7(o) + d 



(35) 



3.2.2 Update of g (e) in multi-bit CS 

In multi-bit CS we have 



q{e) ocexp|(lnp(z|a;,e))g(^)|p(e) 



1 -2/|i „ 4„||2 



«exp|--fT ^||z-e-Aa3||2) We(2?e) (35) 

oc exp |-^a-2 lie - (z - A//)||2| /^ (Pe) 

where /e (Pe) is an indicator function that equals to 1 if e G 2?e or otherwise. Hence, q (e) is the product of 
PDFs of truncated Gaussian distributions, i.e., for each m = 1, • • • , M, q {em) is the PDF of a truncated Gaussian 
distribution. As a result, the first moment of Cm, m = 1, • • • , M, used in ( [32] ) can be given in closed form after some 
derivations using the PDF </) (•) and cumulative distribution function (CDF) $ (•) of a standard Gaussian distribution: 



(t>{le^) -4>{Ue 
'^{UeJ-'^ik 



(em)q(e) = (^ T^rFZ. — ^ T^TTT^ + ^em : (37) 



where /Xe„ = [z - Afi)^, U^ and -Ue^ satisfy that Pe,„ = Wh^ + Me,„,o-iie„ + ^iej\ with Ve^ denoting the 
domain of em, {u) = ^= exp ] — ^ [ and <I> {u) = J"^ (j) {t) dt for n € M. 

Remark 4 Consider the case where q (e) is given the complete flexibility in 1-bit CS. Note that entries of a point 
in T)(, are no longer independent of each other in such a case, leading to that q (e) is the PDF of a truncated multi- 
variable Gaussian distribution with e constrained in a nonconvex set D^ defined in ([9]). As a result, the calculation 
of (e) /g\ is in general computationally intractable in our considered CS problems where the dimension ofe is large. 



3.2.3 Update of q (e) in 1-bit CS 



According to |[8|, this is equivalent to finding an MAP estimator of e with its posterior distribution defined in (36 1. 
So we have 

(e),(e) = e° (38) 

with 



e° = ar, 



^e™RM i ^^^ i ~2^ ^ "^ - (^ - AfJ')\\l \ le {Ve) \ 



= arg min ||e — [z — A^)\\2 

eeVe 

= Vv. {z - A^l) 

where V^ is defined in (J9]l and Vv (f) denotes a projection of a point v onto a set V. The calculation of Vve. (") with 
the nonconvex set Pg is provided in the following lemma. 

Lemma 1 For a vector v G M , let v = —sgn (z) v. Denote X the index set of all positive entries of v. Let 
X^ be its complementary set. IfX is nonempty, then let e* G M with ej = nj^ and ejc = 0. Otherwise, let 

e| = < *° ' i, • ' /"'' ^ = 1) ■ ■ ■ ) M, where io € {!,■ ■ ■ , M} with Vi^ = max {v). Then we have 
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0, otherwise, 

e* = Vve i"^) with Df, as defined in J9 




Figure 2: An illustration of Lemnia[T]with nonnegative entries of e. The unit circle in the first quadrant composes 
of Pe- Projections of four possible v's are shown. 

Proof See Appendix. 

Lemma[T]tells how to calculate the projection onto the nonconvex set V^, defined in (|9]l. An illustration of Lemma 
[T]is presented in Fig. [2} where we consider the two dimensional case with both entries of e nonnegative. The unit 
circle in the first quadrant composes of V^. Projections of four possible v's are shown. The resulting algorithm is 
summarized in Algorithm [T] named as variational message passing with quantization (Q-VMP). 



Algorithm 1: Q-VMP 



Input: sensing matrix A, quantized measurement z, domain of quantization error Dg, and noise variance a 



n 



1. initialize (a„ ^) , -, 

\ " / g(oi) 

2. while not converged do 

3. update SI according to (33 1; 
4. 

5. 
6. 

7. 



h--- ,N,{v),Miind{e) 



Iqie)' 



update fx according to (32 1; 
update (a,^^>g(^) and (a„)^(„), n 
update {'n)gir]) according to (35 1; 
update {e)q(e) according to (37 1 



, A^, according to ( 34 1; 



8. end while 

Output: recovered signal x 



I in multi-bit CS and (38 1 in 1-bit CS respectively; 



/x. 



3.3 The Noise Free Case 



In this subsection we consider Q-VMP in the noise free case. We first consider the data consistency. A consistent 
recovery means that the observation can be reproduced from the recovered signal. Empirical results suggest that 
a consistent recovery result in less errors |T6j22|. A theoretical proof is provided in [23] on the 1-bit case. The 
following analysis applies to both multi- and 1-bit CS. Taking o"^ — )• at both sides of ( 32 1 and ( 33 1 gives 



Thus we have 



/^ 



A 



-1 



_ 1 



A(A 



-1^ 



2 
q{a) 



t 



e), 



A 



'q{a) 



A 



-1 



_ 1 



A(A 



-1 



_ 1 



An —fZ — e^V. 



y 



t 



A (A 



'q{a) 
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i.e., Q (A/j,) — )• z, which indicates that the recovered signal reproduces the observation at each iteration. 

We next consider the update of q (e) in such a case. As cr^ —)• we see that q (e) degenerates into a single-point 



distribution by (36 1, that coincides with the stricter assumption in (31 1 in 1-bit CS. 



3.4 Pruning a Basis Function 

The most difficult computation of Q-VMP is the calculation of I] that is the inverse of an A^ x A^ matrix. Using the 
Woodbury matrix identity, we have 

^ = (A-^>i) - {^-xu ^"^-"^ (^-%u 

with C = a^I + A (A~^)~^, A^ being an M x M matrix. Hence, to calculate S needs O (min {iV^, N'^M]) 
operations. It is noted that if Q-VMP produces some (a~^) , , — )• +oo with n G {1, • • • , A^}, then the correspond- 
ing basis An can be removed from the model. To further speed up Q-VMP, we prune a basis An from the model (to 

) . s is larger than a certain threshold Tpmning- 



reduce A^) when the corresponding parameter (a„^) / -.is larger than a certain threshold Tpmning- Similar basis 



pruning approaches have been used in \^^- 

4 Numerical Simulations 

In this section, we study the performance of the proposed observation model and Q-VMP algorithm in comparison 
with existing ones by numerical simulations. 

4.1 Experimental Setup 

Model Parameters in Q-VMP: We set e = 0, c = 1 and d = Oin the Gaussian-Gamma-Gamma prior as discussed 
in Subsection l2.2l 

Quantizer: In multi-bit CS, a uniform unsaturated quantizer is defined in (Q with L = 2^ , equispaced 

uo,ui- ■ ■ , ul and Vi = {ui + Uj+i) /2, i = 1, • • • , L — 1. In addition, we let u^ = \\y\\^ and uq = — \\y\\^ 
in each trial. For a saturated quantizer, we set uq = — c«, ul = +oo. 

CS problem generation: In our experiment, we set A^ = 500, K = 10, and vary the bit budget (total bits of 
all measurements) in {50, 100, • • • , 1000}. In each trial, a K-sparse signal of length A^ is generated with Gaussian 
distributed nonzero entries and then scaled to unit norm. Entries of the sensing matrix A are generated independently 
according to a Gaussian distribution J\f (O, M~^). Thus the noise free measurement y^ = Ax has unit norm in 
expectation. To obtain a desired SNR, a white Gaussian measurement noise n is added with the noise variance 
cj^ = Af~^10~Tfr. Then the quantized measurement z = Q (y) is preserved for the following signal recovery. 

Performance metrics: Three metrics are considered, including reconstruction SNR (RSNR), sparsity level of 
the recovered signal and computational speed. RSNR is defined as 

RSNR = -201ogio \\x - SII2 

where x denotes the recovered signal of x. The sparsity level is measured by the support size of the recovered signal. 
The computational speed is measured by the CPU time usage. All results are averaged over 200 trials. 

4.2 Model Efficiency 



We first study the efficiency of the observation model in ( 10 1 introduced in this paper for quantized CS. We consider 



the multi-bit CS problem with a uniform unsaturated quantizer as an example. In existing methods that account for 



measurement noise, e.g., in 1 17 1, the quantization error and the noise are typically coupled and treated as a Gaussian 
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Figure 3: Reconstruction SNRs of VMP algorithms implemented respectively based on the proposed observation 
model in ([TOJI, an existing one that couples the quantization error and measurement noise, and conventional CS 
(oracle-aided quantized CS) as an upper boundary. 



noise (only the energy information is used). Then the quantized CS problem is transformed into a conventional 
one. We refer to this formulation as existing method hereafter. In this subsection we compare the signal recovery 
performance of the proposed formulation in ([TO]) with the existing one. Naturally, we use the proposed Q-VMP 
algorithm with our formulation. A corresponding algorithm for the existing formulation is thus VMP introduced 
in |[9| for conventional CS. The latter algorithm can be considered as a simplified version of Q-VMP with the 
quantization error e fixed throughout the algorithm. In addition, we consider the performance of conventional CS 
for comparison where the true-valued measurements are used. The conventional CS problem can be considered 
as one in quantized CS by using oracle information of the quantization error. So its performance acts as an upper 
boundary of the quantized CS problem. 

In our experiment, we set SNR = 30dB and the bit depth B = A that leads to the number of quantized measure- 
ments varying from 12 to 250. Q-VMP is implemented as follows. 






0. We set r, 



pruning 



where a 



vJ-l| 



-1\-1 



1/ \Alz\, n = 1, • • • , iV, (r/)^(^) = 1 and (e)^(^) 

— < 10~^ or the maximum number of iterations, set to 2000, is reached, 

T 



Q-VMP: We initialize (a"^) 
10^. Q-VMP is terminated if 

{oii^} , ),''' j{'^N^} f ) and the superscript J indicates the iteration. 
The VMP algorithm for the other two cases is similarly implemented. The true noise variance is used in Q- 
VMP and conventional CS. For VMP with the existing formulation, we set the noise variance to r^/12 -|- o"^ where 
r = ui — uq denotes the quantization bin width. This noise variance corresponds to a Gaussian noise whose 
energy is comparable with that of e -|- n under the assumption that e is uniformly distributed and independent of n. 
Reconstruction SNRs of the three methods are depicted in Fig. |3] It can be seen that the VMP algorithm based on 
the proposed observation model is consistently better than that with the existing formulation though there is still a 
large gap between it and the oracle-aided case. 
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Figure 4: Performance comparison of Q-VMP, BPDN and LIRML with bit depth B = A. SNR 
Averaged reconstruction SNR; (b) Averaged support size of recovered signal; (c) Averaged CPU time. 



30dB. (a) 



4.3 Performance Comparison in Multi-bit CS 

4.3.1 Unsaturated quantizer 

In multi-bit CS, we first consider the case of a uniform unsaturated quantizer. As in the last subsection, we set 
SNR = 30dB and 3 = 4. Besides Q-VMR we also use BPDN (3) and LIRML |[18| to recover the signal for 
comparison. Q-VMP is implemented as in the last subsection. The other two are implemented as follows. 

BPDN: BPDN solves the problem in (|2]) with y replaced by z and is implemented using £i-magic 1 33 1. We let 



e = \\z — Ax\\2 in our implementation to achieve the best result though it is unavailable in practice. 
LIRML: LIRML solves the problem 

minjAplli -log /ml (Ax)} 

X 

where /ml (•) is the likelihood function of the observation with — log /ml (^^) being a convex function of x, and 
A > is a regularization parameter. In general, a larger A leads to a sparser solution. Since there are no available 
guidelines for the choice of A so far, we choose A such that the recovered signal has the optimal RSNR. Additionally, 
we set r = ..'^.,,2 , e = 10^^ and /3 = 0.5. Readers are referred to [181 for their interpretations. 

^ 1—1 i , j 

The experimental results are shown in Fig. ^ where red solia lines denote Q-VMP, black dashed dot lines 



denote BPDN, and blue dashed lines denote LIRML. Fig. 4(a) depicts the averaged reconstruction SNRs of the 
three algorithms. A significant improvement of the reconstruction SNR can be observed using the proposed Q-VMP. 



It is over 6dB in comparison with LIRML and about an amplitude for BPDN. Moreover, Fig. 4(b) shows that Q- 
VMP produces the sparsest solution. It is noted that LIRML can produce a sparser solution by setting a larger value 



of A as in 1 18 1 but at the cost of a lower RSNR since it produces the optimal RSNR in our implementation. Fig. 4(c) 
shows that the speed of Q-VMP is comparable with that of BPDN and LIRML. Implemented with the basis pruning 
approach, Q-VMP is faster when more measurements aie acquired since it is observed in such a case that the basis 
pruning approach works more efficiently. 



4.3.2 Saturated Quantizer 

We next consider the case of a saturated quantizer. We adopt the same experimental setup but a saturated quantizer 
where a noisy measurement falls in each quantization interval with the same probability. Since both the sens- 
ing matrix and measurement noise are Gaussian in the experiment, the noisy measurements are i.i.d. Gaussian 
Af (0, M~^ + CT^). Then it is easy to get the quantizer. As a result, 12.5% of the measurements are saturated in 
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Figure 5: Reconstruction SNRs of Q-VMP and LIRML with a saturated quantizer, as well as those with the unsatu- 
rated quantizer in Fig. [4] 
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Figure 6: Performance comparison of Q-VMP, BIHT and CVXP in 1-bit CS. SNR = lOdB. (a) Averaged recon- 
struction SNR; (b) Averaged support size of recovered signal; (c) Averaged CPU time. 

expectation. BPDN is inappropriate in such a case. We compare Q-VMP only with LIRML. The averaged recon- 
struction SNRs of Q-VMP and LIRML are presented in Fig. |5](red solid lines). Q-VMP obtains a RSNR of about 
lOdB higher than LIRML when sufficient measurements are acquired. The performance of the two algorithms on 
support size and speed is similar to that in the uniform quantizer case and is omitted. 

The experiment above may shed light on the optimal quantizer design for Q-VMP. By comparing the perfor- 
mance of Q-VMP in the two quantizer scenarios, it can be seen from Fig. |5]that the saturated quantizer outperforms 
the uniform unsaturated one when more measurements are taken for Q-VMP while it is not so clear for LIRML. We 
pose the problem of the optimal quantizer design for Q-VMP as a future work. 



4.4 Performance Comparison in 1-bit CS 

The bit-depth B = \'\rv 1-bit CS. We set SNR = lOdB. In such a case, 9.75% measurements flip their signs due 
to the noise in expectation. We compare Q-VMP with the state-of-the-art algorithms BIHT |[23| and the convex 



programming approach in p5| , denoted by CVXP. The three algorithms are implemented as follows. 



Q-VMP: We initialize (a"^) 



q(pL) 



Mj 



,sgn(zj 



,n 



l,---,iV,(r/) 



iWi 



1 and (e) 



q(e) 



-sgn (z) jyfM. 



15 



As addressed in Remark[3j the effective noise level in 1-bit CS is much lower than the true one. We empirically find 
that it is a good choice to set the noise variance in Q-VMP to lO^'^o-^ (Q-VMP is slow in the case of a very small 
noise variance, which is the reason why we consider a lower SNR in 1-bit CS). We set Tpmning = 10^ and terminate 
Q-VMP as in multi-bit CS. The recovered signal is finally scaled to unit norm for comparison with the original one. 
BIHT: The oracle information of K is used in the hard thresholding operation, i.e., BIHT is certain to return a 
reconstruction with K nonzero entries. BIHT is terminated if the Hamming error of the current recovery is below 
the expected Hamming error or the maximum number of iterations, set to 1000, is reached. Readers are referred 
to ||23[ for the definition of the Hamming error. 



CVXP: CVXP solves the convex problem in (26 1. It does not require the noise information but needs to know 



the signal spai^sity level. We set 5 = ||a;||^ in (26 1 (using the oracle information of ||a;||j^) to achieve the best result. 



We implement it using CVX |34| (its speed can be accelerated using faster solvers). 



Our experimental results are presented in Fig. ml In all figures, red solid lines denote Q-VMP, black dashed 



dot lines denote BIHT, and blue dashed lines denote CVXP. It is shown in Fig. 6(a) that the proposed Q-VMP 
outperforms consistently the other two algorithms in the recovery accuarcy. From Fig. 6(b) it can be seen that Q- 
VMP produces a sparser solution than CVXP on average while BIHT uses this oracle information. Fig. 6(c) shows 
that the computational speed is an disadvantage of Q-VMP. 

5 Conclusion and Future Work 

The sparse signal recovery problem from noisy quantized compressive measurements was studied in this paper. A 
unified framework was presented that applies to both multi- and 1-bit CS. Moreover, the unified framework incor- 
porates the unsaturated and saturated cases in the multi-bit quantization. Relations between the proposed unified 
framework and existing quantized CS methods were studied and a unified algorithm for quantized CS was pro- 
posed based on variational Bayesian inference. Numerical simulations were carried out, showing that the proposed 
framework and algorithm improve the signal recovery accuracy in comparison with the existing results. 

A convex formulation of the noisy 1-bit CS problem has been studied in f^S] with guaranteed signal recovery 



performance. This paper has introduced a different convex formulation (problem (18 1 with f (x) = \\x\\-^ and 
s = 1) that explicitly exploits the noise information and does not need the knowledge of the signal sparsity. One 
future work is to explore its theoretical guarantee. One drawback of Q-VMP is its high computational complexity 
due to an inversion of a high dimensional matrix at each iteration though it has been greatly alleviated with a basis 
pruning approach adopted in this paper. Q-VMP may suffer from computational issues in the case of very high 
dimensional problems. Thus another future work is to develop fast alternatives to the current implementation. Since 
the signal recovery accuracy in multi-bit CS is very different when a different quantizer is adopted, as shown in 
the present paper and in [17], to design the optimal quantizer that minimizes the signal recovery error is another 
interesting future research topic. 

Appendix 

Proof of Lemma [I] 

To prove Lemma [T] we first show the following three lemmas. 

Lemma 2 Assume w G M™, w )^ and w ^ 0. Then v* = ,,^,, w is an optimal solution to the problem 

max f iv) = w^v, subject to I ^^ ' (39) 

with scalar C > 0. Consequently, f* = f {v*) = C ||t«||2. 

16 



Proof By the Cauchy inequality we have / (v) < \\w\\2 \\v\\2 = C ||w||2- The equality holds if v is in the form of 
V*. So, we get the conclusion. 

Lemma 3 Assume w G R"^, w ^ 0. Let v* € M"* with «* = <„, °' far i = I,--- ,m, where 

[0, otherwise, 

io £ {1, • ■ ■ ) "-} with Wig = inin (w). Then v* is an optimal solution to the problem 

T 1 ■ f Il'l'll9 = C, 



mill f (v\ =w V, subject to < „ ' (40) 

with scalar C > 0. Consequently, f* = f (v*) = Cmin (w). 

Proof The case of m = 1 is obvious. We next prove the lemma in the case of ?n, = 2 and then use induction to 
complete the proof. Suppose m = 2. We need to show that if wi < W2, then v^ = C and V2 = 0. By substituting 

vi = a/C^ — V2 into / (i?), we have 



g{v2):=f(^^/c^^lv2^=^/' 



C^ — V2W1+ V2W2- (41) 



It is easy to show that g' (V2) > if < t;2 < ifnr > and g' (V2) < if pf- <V2<C. So, the minimum of g M 
can only be achieved at the boundary of the interval [0, C]. By g (0) = Cwi, g (C) = Cw2 and wi < W2, we have 
v*=0,vl = C and /* = Cmin (w). 

Suppose the lemma holds in the case m = n — 1 with n > 3. We next show that it holds in the case m = n. 
Denote i?_i = [v2, ■ ■ ■ ,Vn] , if-i = [w2, ■ ■ ■ ,Wn] ■ Then we have w'^^v^i > ||i7_i||2 min (tf;_i). By ||t'-i||2 = 
y^C^ — vf we have 

/ (v) = wivi + i(;-^]^t;_i 

> wiVi + min (tu-i) v C^ — f 1 (42) 

> Cmin (w) . 

It is obvious that the equality holds if v is in the form of v*. 

Lemma 4 For w G R"*, denote I the index set of all positive entries ofw. Let X^ be its complementary set. IfX is 

nonempty, then let v* G M™ with v%- = ,, '^\ and vX-c = 0. Otherwise, letv* = <,'' , . ' for i = 1,- ■ ■ ,m, 
^ ^ W'^xh ^ * [ 0, otherwise, 

where iq ^ {!,■ ■ ■ , m} with Wig = max (w). Then v* is an optimal solution to the problem 

max f (v) = w V, subject to < „ ' (43) 

vm"^ ^ ^ -^ \ 17 :>= 0. ^ ^ 

2 ||2 ii2 I I I I 

Proof Consider first the case where X is nonempty. By ||i'x||2 + ||i'z=|l2 = ll^lb ~ ^ and Lemmas p^ and p^ we have 

< \\vx\\2\\wx\\2- ||t'xc||2min(-wxc) (44) 

< Ikx|l2- 

It is obvious that the equality holds if v is in the form of v*. 

In the other case where I is empty, we have w ^ 0. By ||i>||2 = 1 and Lemma [3] we have 

/ (v) = — (— lu ) V < — \\v\\2 min {—w) = max (w) . (45) 

The equality holds again if v is in the form of v* . 
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Proof of Lemma [T| To show e* = Vv^ ('^)» i-^-' 



e = are; mm e — v\ 



we need only to show 



2 ' 



or 



* T 

e = are max i> e, 



i; e > i; e, for any e G Pe, 



or further 

if (sgn (e) e*) > ^^ (sgn (e) e) , for any e G Pe- 
lt is equivalent to showing that —sgn (z) e* is an optimal solution to the problem 

—T— , • f Hello = 1, 

maxt; e, subiect to < „ 

e ' -" \e:^0. 

This is a direct result by applying Lemma|4] So, we complete the proof. 
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